.ipynb .pdf

repository open issue

Contents

# # Cardiac mechanics Benchmark - Problem 2
# This code implements problem 2 in the Cardiac Mechanic Benchmark paper
#
# > Land S, Gurev V, Arens S, Augustin CM, Baron L, Blake R, Bradley C, Castro S,
# Crozier A, Favino M, Fastl TE. Verification of cardiac mechanics software:
# benchmark problems and solutions for testing active and passive material
# behaviour. Proc. R. Soc. A. 2015 Dec 8;471(2184):20150641.
#
from pathlib import Path

import dolfin

try:
    from dolfin_adjoint import Constant, DirichletBC, Mesh, interpolate
except ImportError:
    from dolfin import DirichletBC, Constant, interpolate, Mesh

import pulse
from fenics_plotly import plot

here = Path(__file__).absolute().parent

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-5-ed67a6c2e314> in <module>
----> 1 here = Path(__file__).absolute().parent

NameError: name '__file__' is not defined

geometry = pulse.HeartGeometry.from_file(pulse.mesh_paths["benchmark"])
# geometry = pulse.geometries.benchmark_ellipsoid_geometry()

2021-11-23 07:27:26,725 [852] INFO     pulse.geometry_utils: 
Load mesh from h5

# Create the material
material_parameters = pulse.Guccione.default_parameters()
material_parameters["CC"] = 10.0
material_parameters["bf"] = 1.0
material_parameters["bfs"] = 1.0
material_parameters["bt"] = 1.0

material = pulse.Guccione(parameters=material_parameters)

# Define Dirichlet boundary. Fix the base_spring
def dirichlet_bc(W):
    V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
    return DirichletBC(
        V,
        Constant((0.0, 0.0, 0.0)),
        geometry.ffun,
        geometry.markers["BASE"][0],
    )

# Traction at the bottom of the beam
lvp = Constant(0.0)
neumann_bc = pulse.NeumannBC(traction=lvp, marker=geometry.markers["ENDO"][0])

# Collect Boundary Conditions
bcs = pulse.BoundaryConditions(dirichlet=(dirichlet_bc,), neumann=(neumann_bc,))

# Create problem
problem = pulse.MechanicsProblem(geometry, material, bcs)

# Solve problem
pulse.iterate.iterate(problem, lvp, 10.0, initial_number_of_steps=200)

2021-11-23 07:27:26,967 [852] INFO     pulse.iterate: Iterating to target control (f_28)...

2021-11-23 07:27:26,968 [852] INFO     pulse.iterate: Current control: f_28 = 0.000

2021-11-23 07:27:26,969 [852] INFO     pulse.iterate: Target: 10.000

([Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 44),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 84),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 111),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 146),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 181),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 267),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 310),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 353),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 396),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 439),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 482),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 525),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 568),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 611),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 638)],
 [Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 42),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 82),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 109),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 144),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 179),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 265),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 308),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 351),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 394),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 437),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 480),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 523),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 566),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 609),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 636)])

# Get displacement and hydrostatic pressure
u, p = problem.state.split(deepcopy=True)

endo_apex_marker = geometry.markers["APEX_ENDO"][0]
endo_apex_idx = geometry.vfun.array().tolist().index(endo_apex_marker)
endo_apex = geometry.mesh.coordinates()[endo_apex_idx, :]
endo_apex_pos = endo_apex + u(endo_apex)

print(
    ("\n\nGet longitudinal position of endocardial apex: {:.4f} mm" "").format(
        endo_apex_pos[0],
    ),
)

Get longitudinal position of endocardial apex: 30.7991 mm

epi_apex_marker = geometry.markers["APEX_EPI"][0]
epi_apex_idx = geometry.vfun.array().tolist().index(epi_apex_marker)
epi_apex = geometry.mesh.coordinates()[epi_apex_idx, :]
epi_apex_pos = epi_apex + u(epi_apex)

print(
    ("\n\nGet longitudinal position of epicardial apex: {:.4f} mm" "").format(
        epi_apex_pos[0],
    ),
)

Get longitudinal position of epicardial apex: 32.0417 mm

V = dolfin.VectorFunctionSpace(geometry.mesh, "CG", 1)
u_int = interpolate(u, V)
mesh = Mesh(geometry.mesh)
dolfin.ALE.move(mesh, u_int)

fig = plot(geometry.mesh, color="red", show=False)
fig.add_plot(plot(mesh, opacity=0.3, show=False))
fig.show()

previous

<no title>

next

<no title>